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We explore the available parameter space that gives reasonable fits to the total charm cross section 
to make a better estimate of its true uncertainty. We study the effect of the parameter choices on 
the energy dependence of the J/tp cross section. 

I. INTRODUCTION 

Because the charm quark mass is finite, the total charm production cross section can be calculated in perturbative 
^ ■ QCD. However, the charm quark mass is relatively light so that there are large uncertainties due to the choice of quark 
^ , mass, factorization scale and renormalization scale Typical lovi^er limits of the factorization and renormalization 
■ scales are half the chosen charm quark mass [l]- Here, the parton densities are subject to backward evolution since 
the factorization scale is below the minimum scale of the parton densities. In addition, for renormalization scales 
below 1 GeV, the strong coupling constant as becomes large and the perturbative expansion is unlikely to converge. 
Thus is it worth evaluating our assumptions concerning these parameters to determine whether we can find a set of 
physically defensible mass and scale parameters that reduce the cross section uncertainty. 

In the color evaporation model of J/ ij) production, the J/ ip cross sections are calculated with the same set of mass 
and scale parameters as open charm production Q . As we will show, the parameters used to calculate the uncertainty 
on the charm quark cross sections in Refs. 0,01 do not place stringent bounds on the J/ip production cross section. 
We therefore seek to place limits on the J/-0 cross section calculated in the color evaporation model for the first time. 

In this paper we will explore the charm quark mass and scale parameter space to reduce the uncertainty on the 
charm total cross section. In Sec. [Hi we use existing data to set limits on the factorization and renormalization scales. 
We will then calculate the lepton distributions from heavy flavor decays with our parameter limits and compare to 
data from RHIC and the LHC as a reality check. Section IIIII describes how we use our open charm results in color 
evaporation model calculations of J /ip production to determine the uncertainties on quarkonium production, both as 
a function of incident energy and as a function of the kinematic variables. In Sec. IIVI we summarize our results. 
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II. SETTING LIMITS ON THE TOTAL CHARM CROSS SECTION 



In our previous efforts to place uncertainties on the total charm cross section using the same fiducial parameter set 
as the FONLL calculation, we found a wide uncertainty band that grew larger at high center of mass energies P, 01 ■ 
At high energies, the lower limit was determined by a factorization scale of half the central value of the charm quark 
mass, assumed to be 1.5 GeV. In this region, the behavior of the gluon density at low momentum fraction, x, when 
the parton densities are backwards evolved caused the cross section to grow unphysically slowly with energy. The 
upper limit of the uncertainty is obtained when the renormalization scale is equal to half the charm quark mass. Here 
the two-loop evaluation of the strong coupling constant gives as > 0.5, too large for convergence of the perturbative 
expansion. Here we will use the total charm cross section data to obtain a more physically motivated set of parameters 
for charm production. 

We calculate the total hadronic charm production cross section in a pp collision directly using the next-to-leading 
order (NLO) matrix elements @ for the total partonic cross section, a, 

aAB{Vs-m'^)= ^ Jdxidx2ff{xi,fi%)f^{x2,fJ.p)aij{s,m^,fi%,fi%), (1) 

where m is the charm quark mass, the factorization scale, is the renormalization scale, -y/s is the partonic 

center of mass energy, xi and X2 are the fractions of the parent proton's momenta carried by the colliding partons 
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and /f are the proton parton densities. The NLO calculation remains the state of the art for the total cross section; 
there is still no complete NNLO evaluation of the total cross section, especially at energies where -y/s ^ m. We use 
the central CTIO parton density set but will also show the variation in the cross section based on all 52 variants 
of the Hessian uncertainty matrix. 

Since Eq. ([Ij is independent of the heavy quark kinematics, it is typical to take hb.,f — m as the central value 
and vary the two scales independently within a 'fiducial' region defined by fifj^F/m with 0.5 < ^iR^p/fn < 2 and 
0-5 < IJ'rI I^F < 2. In earlier work, we used the following seven sets: {{fiF/m, fin/m)} = {(1,1), (2,2), (0.5,0.5), 
(1,0.5), (2,1), (0.5,1), (1,2)} [HIl]. The uncertainties from the mass variation and the combined scale variations listed 
above were then added in quadrature. The envelope containing the resulting curves. 



^ccnt)'^ ~t~ (f^7n,max ^ccnt)^ ; (2) 



^inin '-^ccnt ^ {^fi^min ^ccnt) ~t~ ((^7n,min f^ccnt) j (3) 

defines the uncertainty on the total cross section as a function of center of mass energy. Here (iccnt is the cross section 
calculated with the central set, {fxp /m, fij^/m) = (1, 1) and m = 1.5 GeV, while CTi.max and CTi^min are the maximum 
and minimum values of the cross section for a given mass (i — m) or {^p/m, fxn/m) set in the fiducial region (i = fi). 
Although Eqs. ([2]) and ([3]) have been written for the total cross section, the corresponding maximum and minimum 
values of the differential distributions can be written similarly [2] . 

The charm quark mass we employ in our calculations is the Particle Data Group (PDG) value based on lattice 
determinations of the charm quark mass in the MS scheme at /x = to: m{m) = 1.27 ± 0.09 GeV 0. The fiducial 
cc parameter sets used in FONLL calculations employ a higher charm quark mass, m = 1.5 GeV. None of these 
fiducial parameter sets give a particularly good fit to the total charm data. When Uf = 3 flavors are used, as is proper 
for charm production, there is a wide uncertainty band on CTtot, especially at the center-of-mass energies appropriate 
for colliders, > 200 GeV, primarily due to unconstrained gluon densities at low x for /jp/m = 0.5 < fio/m where 
fiQ is the minimum scale of the parton densities. Previous calculations with lower charm quark masses but higher 
scales [3] agree better with data while avoiding backward evolution of the gluon density at low x. This bias of lower 
masses with higher scales allows us to reduce the uncertainty in the charm production cross section. 

In principle, fitting the data is somewhat problematic since we neglect unknown next-order uncertainties. This is 
particularly true for charm where the mass is relatively small and 0{a'l) corrections could be large. Indeed, approx- 
imate NNLO calculations show that, while the scale dependence is reduced, the K factor between the approximate 
NNLO and the NLO results is similar to that between the NLO and LO calculations [§, [l3]- Since a full NNLO 
calculation is not yet available, we feel a fit that narrows the uncertainties at collider energies is useful, keeping in 
mind that a full NNLO calculation might yield a good fit to the data with higher masses and somewhat lower scales. 

For a fixed charm quark mass, we fit the factorization and renormalization scale parameters to a subset of the 
total charm production data. We use part of the fixed-target data measured with incident protons at beam energies 
Ebci>,m = 250 [III, 360 [13 , 400 [13 , 450 [i3 , 800 [il,[i3, and 920 GeV [l^. We do not include incident pion data in 
the analysis because there have been no new global analyses of the pion parton densities since 1999 ^Ts^ and none of 
the past pion fits are compatible with modern proton parton densities. The Lexan bubble chamber (LEBC) was used 
in the measurements of the NA16 [l3, NA27 [13 and E743 [13 Collaborations. LEBC allowed direct observation of 
the charm production and decay vertices. The first two measurements were made at CERN [l^, [13 while the last 
was made at Fermilab [l3. The 800 GeV E653 measurement at Fermilab used an emulsion target to measure the 
primary production vertex and at least one decay vertex contained within the emulsion volume [l6[. While none of 
these experiments had very high statistics, their results were very clean. The E769 Collaboration used silicon vertex 

detectors to reconstruct D meson (D^, D°/D and Df) decays [Tlj. The NA50 data at i^beam = 450 GeV were 
obtained by studying the lepton pair invariant mass continuum over a range of nuclear targets. The continuum was 
assumed to be a superposition of dimuons from the Drell-Yan process and semileptonic decays of open charm. Since 
the A dependence of open charm and Drell-Yan production is compatible with a linear growth, cFpA_= o'ppA^ the 
charm cross section was obtained from a global fit to the four targets studied (Al, Cu, Ag and W) 'iJ]. The data 
from Refs. 

[iMl[il[i3 

were evaluated in the review of Lourengo and Wohri and adjusted to the values we employ 
in our fits using the most up-to-date branching ratios for the measured decay channels [l9j . 

We also include total cross section data at y/s = 200 GeV from RHIC. There are data from both PHENIX 
and STAR |^21r,23] . The PHENIX measurement is based on inclusive single electron pT distributions in pp collisions 
in the pseudorapidity interval I77I < 0.35. The 'non-photonic' electrons, assumed to come from heavy flavor decays, 
were extracted from the total electron spectrum by subtracting 'photonic' (background) sources. The shape of the 
resulting px distribution is described by a superposition of charm and bottom contributions. The charm contribution 
was extrapolated to pt = to obtain the total charm cross section, O.bblto^ll mb [gCj]. The first STAR data point 
was extracted from d+Au collisions by two independent measurements [2lj. They directly reconstructed Z?'^ — >■ K~^Tr~ 
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Fitted Data 


lip /m 


flR/m 


xVdof 


fixed-target only 

+ PHENIX 
+ STAR (2004) 
+ STAR (2011) 
+ STAR (2012) 


1 1+1.00 
^•^-0.40 

J^-O-0.56 
O c+2.73 
^•0-1.35 
9 1 +2.55 
^•^-0.85 
9 1 +2.21 
^- J^-0.79 


^•"-0.08 
1 c+0.09 
^•"-0.13 

^•"-0.10 

1 6+°" 

^•"-0.12 

1 6,+°'^° 

J^-O-0.11 


1.03 
1.03 
1.53 
1.16 
1.06 



TABLE I: The factorization, fiF/m, and renormalization, fj,R/m, scale uncertainties obtained by fitting subsets of tfie total 
charm cross section data with m = 1.27 GeV. 

decays with < 1 and 0.1 < pr < 3 GeV. STAR also used inclusive non-photonic electrons to study semileptonic 
decays of charm. The initial result, 1.4 ± 0.2 ± 0.4 mb [21], was significantly higher than the PHENIX result but 
compatible within systematic uncertainties. After a reanalysis of the non-photonic electron data and new D meson 
measurements, the STAR cross section reported at Quark Matter 2011 came down to 0.949 ± 0.365 mb [l^] while the 
final result, obtained after our analysis was finished, is reduced to 0.797 ± 0.210 j:[5:262 

mb [2^. This does not change 

the central value of the scale parameters, only reduces the one standard deviation limits. While the final STAR result 
is still higher than the PHENIX cross section, the two results are now comparable within uncertainties. 

We have made five different fits to combinations of the data just described: the fixed-t arg et data (TT| - [l3l ITBl - flTj 
only; adding only the PHENIX data [20*1; adding the PHENIX and the 2004 STAR result [21] ; and finally, including 
the 2011 STAR result [22] and, subsequently, checking how much the results changed when the final 2012 STAR point 
[2^ was added to the PHENIX data. The experimental uncertainties used in the fitting were obtained by adding the 
statistical and systematic uncertainties in quadrature. 

In our analysis, the total charm cross sections were calculated for a range of charm quark masses between 1.18 and 
1.54 GeV in steps of 0.03 GeV. At each mass, we varied fxp/'m between 0.45 and 10.65 while simultaneously varying 
fiR/m between 0.5 and 2.9. The step size in fj,p/m and ^r/tti was 0.05 in both cases. The x^/dof for each parameter 
set was evaluated by comparing the calculated cross sections with each of the five subsets of the data considered. 

The best fit values of /m and ^R/m are rather sensitive to the charm mass. In general increasing the quark mass 
above 1.27 GeV decreases both fip/m and fiR/m. It also tends to increase the x^/dof for each fit. If one plots x^/dof 
for a given value of m as a function of either ^p/m or pLR/m while the other scale parameter is held fixed, typical 
parabolic shapes with a minimum are found. The parabolas grow narrower as the mass increases. When to < 1.2 
GeV, single variable parabolas of x^/dof are rather broad and prefer high fiR^p/m values. For to > 1.5 GeV, the fits 
give ^_f._r/to < 1, close to the minimum ^p of the parton densities and in a region where as(/i|j) is rather large. 

Because the charm quark mass was assigned the value of to = 1.27 ± 0.09 GeV by the PDG, we decided to add 
a penalty to the x^ equal to (to — topdg)^/^Wpdg- With this penalty for deviations from the PDG value of the 
charm quark mass, the minimum x^/dof when varying the charm quark mass within our chosen range was found for 
TO = 1.27 GeV, the PDG mass, for all five subsets of the data. 

The best fit results in all cases are given in Table HI The x^/dof for each fit is also shown. The largest x^/dof is 
obtained when the 2004 STAR point is used since it is high relative to the -^/s dependence of the other measurements. 
We note that the values of fj,p/m found with the later STAR results is more in line with physical arguments than 
that obtained with the 2004 data. 

The uncertainties in the fitted parameters were evaluated from the x^ distributions. We show the fit contours 
in Fig. [1] for the four cases represented in Table ID The x^ contours in ^p/m (a;-axis) and iiR/m (y-axis) are depicted 
at Ax^ — 0.3, 1 and 2.3. The one standard deviation uncertainty in the fitted value of jip/m {^R/m) was taken as 
the maximum extent of the A^^ = 1 contour along the ^p/m {fj,R/m) axis. These uncertainties are included with 
the best fit parameter values in Table HI The one standard deviation uncertainty in the total cross section is the range 
of cross sections resulting from all combinations of fip/m and ^r/iti contained within the A^^ = 2.3 contour. The 
Ax^ = 0.3 contour is included only to guide the eye. 

Using the final STAR data point [23j in the fitting results in the same optimum parameter values for /i p jm and 
fjLR/m. However, the uncertainties on the parameter values are somewhat modified. The upper and lower limits on 
jip/m are reduced by 8% and 4% respectively, while the limits on jJiR/m change by less than 1%, see Table HI Since 
the analysis for this paper was completed before the latest STAR charm data release, we used the limits obtained 
with the preliminary point in our further analysis. 

Note the narrow range in ptRlm relative to the much broader [ip jm range, even for fits to the fixed-target data only. 
Indeed, the largest difference in the fits to the various data sets is in the [ipjm range. The fip/m range compatible 
with the data varies considerably for the different fits, note the difference in ^p jm ranges for the four panels in 
Fig- EI The fixed-target data probe a region of relatively large parton momentum fractions, x ~ 2m/ y/s, equivalent to 
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0.06 < X < 0.12 for 19.4 < -y/s < 40 GeV. This range of x is near the pivot point of the gluon distribution, xg{x, is^), 
as a function of x for a range of factorization scales. The fixed-target data are therefore rather insensitive to the 
evolution of the gluon density as a function of /Xi? so that the results skew toward rather low values of /Xi?/m. 
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d) PHENIX+STAR(2012) 
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FIG. 1: (Color online) The xV^of contours for (a) fixed target data only, (b) including the PHENIX 200 GeV cross section, 
(c) including the STAR 2004 cross section and (d) including the STAR 2012 cross section but excluding the STAR 2004 cross 
section. The best fit values are given for the furthest extent of the A^^ = 1 contours. Note that while the y-ax;is range in the 
same in all four panels, the a;-axis range varies significantly. 

Including RHIC data in the fit introduces greater sensitivity to the low x region although x ^ 0.012 at midrapidity 
is near the high x edge of the low x regime. The PHENIX point, obtained earliest, has the lowest cross section and 
thus requires a lower factorization scale than when either of the two STAR points are included. The STAR cross 
section from 2004, more than a factor of two larger than the PHENIX cross section, requires the largest factorization 
scale of all the fits. Note the high value, fip/ni ^ 10, required to close the Ax^ = 2.3 contour for this fit. The newer 
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STAR measurements |22|,|23|, based on reconstructed D meson decays in pp collisions, rather than on d+Au collisions, 
gives a lower best fit value of ^p/m than the 2004 cross section but still higher than either the fixed-target data only 
or with only the PHENIX measurement at 200 GeV. 

The value of fip/m is strongly dependent on the data sets used in the fits. The uncertainty on np/m is very 
asymmetric with a 100% or greater upper uncertainty relative to the best fit value. The difference between the lower 
limit of the uncertainty on /ip/m and the best fit value is not as large because there is a much greater change in 
xg{x,fijp) at lower factorization scales than when ^ fig. We finally note that the value of ^p/m has the greatest 
effect on the energy dependence of the total charm cross section. 

The best fit value of fiR/m is the same in all cases and the uncertainty is much smaller than for fip/m. These 
uncertainties are also asymmetric but they typically differ by less than 10%, indicating that fiR/m acts to fine tune 
the magnitude of the cross section. Changing fipi/m changes the total cross section by the same factor at all energies 
and does not affect the energy dependence of the cross section. 

Figure [5] shows the resulting energy dependence of the total charm cross section for the four different fits and the 
corresponding uncertainty based on results using the one standard deviation uncertainties on the quark mass and 
scale parameters. If the central, upper and lower limits of jJtpi.p/m are denoted as C, H, and L respectively, then 
the seven sets corresponding to the 'fiducial' region are {{iip/m, fip/m)} = {(C,C), {H,H), {L,L), {C,L), {L,C), 
(C, H), (H, C)}. The upper and lower limits on the PDG value of the charm quark mass are 1.36 and 1.18 GeV. The 
uncertainty band can be obtained for the best fit sets using Eqs. @ and ([3]). The uncertainty bands are shown for two 
cases: the regular fiducial region and including the most extreme cases {fip/m, fin/m) — {H,L) and {L,H). These 
two combinations give the most extreme values of the cross section because the maximum value of fip/ni produces 
the fastest evolution of the parton densities while the minimum value of ^pi/m resits in the largest values of the 
strong coupling constant with (fip/m, fin/m) = {H,L) while the opposite is true for {fj,p /m, iin/m) = {L,H). The 
difference between the outer magenta curves, which include these extremes, and the cyan curves, which do not, is 
very small. Therefore, it is reasonable to neglect the effect of these extremes. 

Note that the fits all result in an asymmetric uncertainty band for ^/s > 100 GeV. This arises because the uncertainty 
in the fits of jip /m is asymmetric, see TablelH with the upper value significantly higher than the lower. As jip increases 
so that ptp ^ /Lto, the evolution of the gluon density with ^p \s reduced for the upper limit of ^p/m. However, the 
closer the lower limit of the fitted iip is to /io, the stronger the factorization scale evolution of the gluon density 
becomes, giving a greater difference between the central value of iip/m and the lower limit than between the central 
value and the upper limit. 

All the fit results shown in Fig. [2] agree equally well with the fixed-target data. However, the fit to the fixed-target 
data alone gives the lowest cross sections at collider energies, ^/s > 200 GeV. The low factorization scale values result 
in a slowing of the growth of the total cross section. The narrowest uncertainty band is obtained when the 2004 STAR 
measurement is used in the fit because it requires the largest factorization scale. Despite this, the top of the calculated 
uncertainty band does not even touch the bottom of the uncertainty on measurement. On the other hand, the most 
recent STAR measurements are compatible with the upper limits of the uncertainty in the fit values. The stronger 
growth in the energy dependence of the total cross section when the RHIC data are included is due to the requirement 
of a larger value of /Xi;-/??! to fit the data, steepening the slope of the energy dependence at large y^. The dot-dashed 
and dashed curves in Fig. (^Jd) were calculated with the preliminary STAR 2011 point |2a. The black dotted curves 
in Fig. [UJd) show the limits on the cross sections calculated using the final STAR value [23] . The difference in the 
calculated upper limits is 0.77% at 200 GeV and 0.70% at 7 TeV while the difference in the lower limits is —3.36% 
at 200 GeV and —12.65% at 7 TeV. There is a smaller difference in the upper limits due to the relatively smaller 
changes in the gluon distributions at low x, high fip, compared to low x, low iip. 

Finally, in Fig. [IJd) we also show the result for a one standard deviation uncertainty in the total cross section 
obtained from the Ax^ = 2.3 contour in Fig. [TJ The resulting band is narrower than the uncertainty band obtained 
from the scale uncertainties in the region of fixed-target data but is compatible with the scale uncertainties at = 200 
GeV. Since it is based on the energies of the data in the fits, it is not extrapolated to either higher or lower energies. 

Lastly, we have added the 2.76 and 7 TeV total cross sections obtained by the ALICE collaboration in pp collisions 
[2^ . These points were not included in our fits. Only the fits where both PHENIX and STAR data are included, giving 
more weight to the RHIC results, have central cross section values close to the LHC data. While both calculations 
lie close to the data, the for the LHC points is 8.67 with the 2004 STAR point and 3.9 with the latest result. 

The individual components of the uncertainty band for the fit including the STAR 2011 data are shown in Fig.[3ja). 
The uncertainty due to the charm quark mass (solid curves) dominates for < 100 GeV where the scale uncertainty 
begins to become comparable. Indeed, the scale variations at fixed-target energies are contained within the curves 
delineating the mass uncertainty. This is very different from the behavior of the fiducial set based on m = 1.5 GeV 
where the scale variation dominates the uncertainty at all ^/s. As the energy increases, the change in a; ~ 2m/ 
due to the mass has a much smaller effect on xg{x,^'^) than the change in the evolution of the gluon density with 
IJ,p. At higher center of mass energies, the curves cluster according to the factorization scale choice. At the top, with 
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FIG. 2: (Color online) The energy dependence of the charm total cross section compared to data for the fits to (a) fixed target 
data only, (b) including the PHENIX 200 GeV cross section, (c) including the STAR 2004 cross section and (d) including the 
STAR 2011 cross section but excluding the STAR 2004 cross section. The best fit values are given for the furthest extent of the 
Ax^ = 1 contours. The central value of the fit in each case is given by the solid red curve while the dashed magenta curves and 
dot-dashed cyan curves show the extent of the corresponding uncertainty bands. The dashed curves outline the most extreme 
limits of the band. In (d), the dotted black curves show the uncertainty bands obtained with the 2012 STAR results while 
the solid blue curves in the range 19.4 < < 200 GeV represent the uncertainty obtained from the extent of the Ax^ ~ 2.3 
contour in the bottom right panel of Fig. [1] 
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FIG. 3: (Color online) (a) The components of the uncertainty band in Fig. [IJd). The central value {m, fip /m, iia/m) = 
(1.27 GeV, 2.10, 1.60) is given by the solid red curve. The solid blue and magenta curves outline the mass uncertainty with 
(1.18 GeV, 2.10, 1.60) and (1.36 GeV, 2.10, 1.60) respectively. The dashed curves outline the lower limits on the scale uncertainty: 
[liF /m, yLR/m) = (2.10, 1.48) blue; (1.25,1.60) magenta; and (1.25,1.48) red. The dotted curves outline the upper limits on the 
scale uncertainty: {fip/m, ^.n/m) = (2.10, 1.71) blue; (4.65,1.60) magenta; and (4.65,1.71) red. The upper and lower dot-dashed 
cyan curves correspond to {^lf /m, ^n/ra) — (4.65, 1.48) and (1.25,1.71) respectively, (b) The uncertainty band on the total 
charm cross section obtained with the FONLL fiducial parameter set centered around (m, fip/m, fiR/m) = (1.5 GeV, 1, 1). The 
central value is given by the solid red curve while the limits of the uncertainty band are shown in the dashed magenta curves. 
The dotted blue curve is the result for (m, /m, iiR/m) — (1.2 GeV, 2, 2). 



the largest growth as a function of -^i are the largest values, ^p/m ^ 4.65. The lowest value of ^_F/m, 1.25, causes 
the slower growth in cross section because the gluon distribution is increasing slowly with decreasing x for this value. 
The uncertainty arising from the range of fiR/m are rather small, due to the narrow range of fit values, and shift the 
overall magnitude of the curves rather than change the slope. 

The spread in the calculations can be compared to the uncertainty band obtained using the fiducial FONLL 
parameter set based on to = 1.5 GeV [H, [|, shown in Fig. ISlJb). The prior by-eye fit to the data using to = 1.2 
GeV, iip/m — iiR/m = 2 0] is also shown in this plot. It gives a better representation of the data than the central 
FONLL parameter set, to — 1.5 GeV, /i^/to = nn/m = 1, and is nearly equivalent to the best x^/dof obtained with 
TO = 1.2 GeV. It also lies rather close to the LHC points (x^ = 23.1) while the central NLO cross section with to = 1.5 
GeV is a factor of ~ 3 below these data (x^ = 142.6). The upper limit of the uncertainty band in this calculation, 
obtained with ^ir/ttl = 0.5, is a factor of ^ 2 larger than the fitted upper limit because as is a factor of ~ 1.6 larger 
for fin = 0.75 GeV than that calculated for the fit results. On the other hand, the lower limit is a factor of 7-8 below 
that calculated with the fit results at LHC energies. At collider energies, the slower growth in the cross section with 
^/s is due to the factorization scale jip — 0.75 GeV {^p I'm = 0.5), below the minimum scale of the PDFs, resulting 
in backward evolution of the gluon distribution. In the fixed-target energy range, the difference is due to the largest 
mass used, m — 1.7 GeV rather than 1.36 GeV. While the uncertainty band obtained with the FONLL fiducial set is 
large enough to encompass all possibilities, it is too wide to give the calculation any predictive power. In addition, 
the scale uncertainty is considerably larger than the mass uncertainty which should not be true for the physical cross 
section. 

The behavior of the gluon density corresponding to the lower, central, and upper values oi for the fits including 
the STAR 2011 cross section are shown in Fig. |4l When jip /ra = 1.25, fip is only ~ 20% higher than so that 
xg{x,^^) is almost independent of x for x < 0.01. As ^p/m increases, the growth of the gluon density at low x 
becomes more pronounced while the uncertainty band becomes narrower for all values of x. It is clear from these 
results that the behavior of xg{x, /ip) determines the growth of the total cross section as a function of center-of-mass 
energy. 

Since the gluon density is not directly measured, the uncertainty in its behavior as a function of x and ^p can be 
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FIG. 4: (Color online) The CTIO gluon distribution, xg{x,^F), is shown for the relevant values of ^p/m for the total cross 
section calculation. The central value of the CTIO gluon distribution is given in the red solid curve while the uncertainty band 
is denoted by the dashed magenta curves. The results are shown for the lower limit of ^p/m, ^p/m — 1.25 (a); the central 
value, np/m = 2.1 (b); and the upper limit, fip/m = 4.65 (c). (d) The corresponding uncertainty on the total charm cross 
section due to the uncertainty in the CTIO gluon distribution is denoted by the dashed magenta lines. The total uncertainty 
due to the mass and scale uncertainty as well as the gluon uncertainty, combined in quadrature, is given by the dot-dashed 
blue curves. 

important. The largest uncertainty can be expected at low scales. To quantify the uncertainty in the gluon density 
that enters into our calculations, we also show the resulting uncertainty band obtained by combining all 52 sets for 
the 26 eigenvectors of the Hessian matrix for the CTIO parton densities. The limits on the behavior of xg{x,^^) 
show the most variation for the lower limit of the factorization scale, see Fig. SJa). There is a sharp increase in the 
upper limit for x < 0.001 while the lower limit on the band has a dip at the same value of x. Using the lower limit of 
jiF/ra — 1.31 from the latest analysis slightly reduces the variation in the band. As iip/ni increases, the growth of 
the gluon density at low x becomes more pronounced while the uncertainty band becomes narrower for all values of 

X. 

The dashed curves in Fig. lU^d) show the uncertainty on the total charm cross section due to the variation of the 
proton parton density. We have used the scale uncertainties from the fit to the 2011 STAR result here. We note 
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that while we have shown the uncertainty bands on the gluon density in Fig. |3Ua)-(c), the cross section uncertainty 
shown here includes the variations in both the quark and gluon densities. In general, the uncertainty due to the 
parton densities is smaller than that due to the scale choice. The combined effect of the mass, scale and parton 
density uncertainties are given by the dot-dashed curves. It is generally only somewhat wider than that due to the 
mass and scale uncertainties alone except for the upper limit of the band at ^/s > 1 TeV. 




5 10 15 0.0 2.5 5.0 7.5 10.0 
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FIG. 5: (Color online) (a) The components of the non-photonic electron spectrum: _B — >■ e (dot-dashed blue); _B — >■ D — >■ e 
(dotted magenta); D e both with the FONLL parameters (solid red) and those for m — 1.27 GeV in Fig. [ijd) (dashed black) 
at \y\ < 0.75 in — 200 GeV pp collisions, (b) The sum of the contributions are compared with the FONLL set for charm 
(solid red) and that with m = 1.27 GeV (dashed blue). The PHENIX data ^ are also shown. 

We now discuss how our results for the mass and scale parameters affect the kinematic distributions of semileptonic 
decays of charm. The state-of-the-art calculational method for single inclusive heavy quark production and decay is 
the fixed-order next-to- leading logarithm approach (FONLL) [2^. In addition to including the full fixed-order NLO 
result , the FONLL calculation also resums [23| large perturbative terms proportional to a" log*^ {pr /m) to all 

orders with next-to-leading logarithmic (NLL) accuracy {i.e. k = n, n — 1). The total cross sections obtained by 
integrating the FONLL kinematic distributions, Eq. should be equivalent to that obtained by convoluting the 
total partonic cross sections with parton densities, Eq. ([T]) when the same number of light flavors are employed. 

The main difference in the two approaches that might affect the total charm cross section is the number of active 
flavors. In the FONLL approach, the heavy quark is treated as an active light flavor at px 3> m. Thus the number of 
light flavors used to calculate as includes the heavy quark, i.e. Tiif + 1 where, for charm, nif — 3 {u, d and s). The same 
number of flavors, nif -I- I, is also used in the fixed-order component of the FONLL calculation for self-consistency. 
Therefore, a total charm cross section calculated in the FONLL approach will automatically be lower than the result 
with the same mass and scale parameters in Eq. m with nn = 3 since as{nif = 4) < as (nif = 3). When the 
renormalization scale is on the order of the quark mass, the difference in the total cross sections at y/s ~ 200 GeV is 
less than 20% However, for iir/itl < 1, asifin) grows faster with decreasing so that the upper limit on NLO 
cross section is up to a factor of two larger than that obtained with FONLL. On the other hand, the lower limit, 
obtained with fia/m = 2, is very similar in the two calculations. Thus whether charm is treated as a heavy (njf) or 
an active (njf -I- 1) flavor in the calculation turns out to be one of the most important influences on the limits of the 
charm uncertainty comparing the NLO and FONLL results. When the total charm cross section is calculated with 
nif in the FONLL approach, i.e. the charm quark is treated as a heavy rather than an active flavor, the results are 
in agreement with the NLO calculations [28]. 

The calculation of the inclusive electron spectrum from heavy flavor decay involves three components: the pr and 
rapidity distributions of the heavy quark Q, calculated in perturbative QCD; fragmentation of the heavy quarks 
into heavy hadrons, Hq, described by phenomenological input extracted from e^e~ data; and the decay of Hq into 
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electrons according to spectra available from other measurements, schematically written as 

E(fa{e) _ EQ(Pa{Q) 



dp-^ dpQ 



■D{Q^HQ)®f{HQ^e) (4) 



where the symbol ® denotes a generic convolution. The fragmentation of quarks into hadrons is denoted D{Q — )• Hq). 
The electron decay spectrum, /{Hq — > e), accounts for the semileptonic branching ratios. 

Figure [5] shows the lepton spectra arising from semileptonic heavy flavor decays at y/s — 200 GeV, all calculated in 
the FONLL approach. The B e and B ^ D ^ e bands, as well as the red — > e band, are calculated with the 
same fiducial set of parameters as in Ref. [5]. The black dashed curves represent the D — )■ e band calculated for our 
best fit parameter set, including both the PHENIX and 2011 STAR total cross sections. The new Z? — e calculation 
is much narrower. It lies completely within the uncertainty band based on the fiducial parameter set with the central 
value of the charm quark mass fixed at m = 1.5 GeV. At low pr, PT < 2.5 GeV, the new D e band is near the top 
of the fiducial FONLL D e band while for pT > 7.5 GeV, the new set gives a result near the bottom of the fiducial 
D ^ e band. The transition from dominance of the electron spectra by charm decays to bottom decays happens at 
lower PT with the m = 1.27 GeV set. 

The right-hand side of Fig. [5] shows the sum of the D e, B e and _B _D e for the two cases. The PHENIX 
non-photonic electron data are compatible with the top of the sum of the uncertainty bands with the fiducial FONLL 
set. However, the data now lie somewhat above the band obtained when the best fit Z? — e contribution to charm 
production replaces the fiducial contribution. The agreement with the data is worst at intermediate values of pt, 
PT ~ 5 GeV, where the c and b contributions to the electron spectra are nearly equal. This discrepancy is not so 
surprising because the best fit parameters were obtained using the NLO QCD calculation with n\{ = 3 fiavors while 
the FONLL calculation, both in the fixed-flavor scheme and in the full FONLL result, uses nif + 1 = 4 since the heavy 
quark is treated as an active light flavor over all pT- The value of asipji) obtained with four light flavors is smaller 
than that obtained with three light flavors, even for the same value of fin. More importantly, the range of factorization 
scales is larger for our best fit case, 1.59 < fip < 5.9 GeV (1.25 < [ipjrti < 4.65) instead of 0.75 < < 3 GeV 
(0.5 < iip/m < 2) for the fiducial set. The higher factorization scales cause the pt distribution to fall off faster with 
Pt in the best fit case. The difference is apparent already in the charm quark pt distributions and would be enhanced 
for the semileptonic charm decays to electrons since the decay leptons carry only ~ 30% of the parent hadron pT [29| 

Figure [6] compares our calculations with the ALICE single muon data in the forward rapidity region, 2.5 < y < 4 
[30| . The data are given for 2 < pT < 12 GeV, both over the full rapidity region. Fig. IH^a), and separated into five 
rapidity bins, each 0.3 units wide, Fig.|6ljb). The calculations with both the fiducial charm parameter set (solid) and 
our charm fit (dashed) are compared to the data in Fig. [Sfa) . The two bands are indistinguishable for pT > 5 GeV. 
Therefore, for clarity, we compare the muon pT distributions in the narrow rapidity bins to only our calculations with 
the mass and scale parameters from the charm fit. The calculations agree well with the measurements over the entire 
Pt range. 

In Fig. in^c) and (d) we present the results as a function of rapidity integrated over the same pt range as the data, 
2 < Pt < 10 GeV. Figure[6jc) shows the upper and lower limits of the FONLL calculations of B ^ ^ and B ^ D ^ jj. 
in the dot-dashed and dotted curves respectively. The FONLL D ^ jj, uncertainty bands with the fiducial charm 
parameter set are shown by the solid curves while the dashed curves are calculated with the charm flt parameters. The 
sum of the heavy flavor decay contributions to the rapidity distribution are compared on a linear scale in Fig. [BUd). 
The pT-integrated ALICE data agree well with both calculations. The results with the fitted charm parameter set 
narrow the uncertainty band without sacriflcing consistency with the measured data. 

While the agreement between the lepton measurements at RHIC and the LHC and our calculations is encouraging, 
as noted here and in Ref. 0, there is significant admixture of semileptons charm and bottom decays, particularly at 
lepton Pt > 4: GeV. A better test of our results would be a comparison to open charm hadron data^ Thus, in Figs. [7] 
and [SI we show the D° (a), D+ (b) and D*+ (c) distributions in the ALICE [32] and the LHCb [33i] acceptances at 
midrapidity and forward rapidity respectively. 

Figure [7] compares the FONLL calculations with the fiducial parameter set (in red) with the fitted parameters based 
on m = 1.27 GeV (in blue). The upper and lower limits of both bands are shown. While the ALICE data are in 
agreement with the upper limits of both calculations, the large D meson uncertainty is reduced at low pT with the 
fitted parameter set. 

Figure [8] shows the upper and lower limits of the FONLL calculation based on m = 1.27 GeV in the five rapidity 
intervals of Ay = 0.5 in the range 2 < y < 4.5 covered by the LHCb detector. In most cases here also the agreement 



" For a complete discussion of LHC predictions using the fiducial FONLL parameter set, see Ref. [Sl 
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FIG. 6: (Color online) Our calculations are compared with the ALICE inclusive single muon data from heavy flavor decays 
[3^ at — 7 TeV. (a) Comparison of the single lepton pr distributions in the rapidity interval 2.5 < y < 4 at ^/s — 7 TeV 
calculated with the FONLL set for charm (solid red) and the fitted set with m — 1.27 GeV (dashed black), (b) The contributions 
to the pt distributions in (a) divided into rapidity bins, from top to bottom: 2.5 < y < 2.8 (solid red); 2.8 < y < 3.1 (solid 
blue); 3.1 < y < 3.4 (dashed red); 3.4 < y < 3.7 (dashed blue); and 3.7 < y < 4 (dot-dashed red). The top curves are shown 
at their calculated value, the others are scaled down by successive factors of 10 to separate them, (c) The components of the 
rapidity distribution at y's = 7 TeV with 2 < pr < 10 GeV, B fi (dot-dashed blue); B ^ D fi (dotted magenta); D ^ fj, 
both with the FONLL parameters (solid red) and those for m — 1.27 GeV in Fig. [ijd) (dashed black), (d) The sum of the 
contributions are compared with the FONLL set for charm (solid red) and that with m — 1.27 GeV (dashed black). 
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FIG. 7: (Color online) Our calculations are compared with the reconstructed ALICE (a) D", (b) Z)+, and (c) D*+ meson data 
[3^ at ^/s = 7 TeV in \y\ < 0.5. The FONLL uncertainty bands with the fiducial charm parameter set are shown by the red 
solid curves while the blue dashed curves are calculated with the charm fit parameters. 
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FIG. 8: (Color online) Our calculations are compared with the reconstructed LHCb (a) D*^, (b) D'^ , and (c) D*"^ meson data 
[3^ at -^5 = 7 TeV in the rapidity intervals: 2 < y < 2.5 (solid red); 2.5 < 3/ < 3 (solid blue); S < y < 3.5 (dashed red); 
3.5 < 2/ < 4 (dashed blue); and 4 < 1/ < 4.5 (dot-dashed red). The curves are calculated with the charm fit parameters. The 
sets of results are separated by a factor of 10 between rapidity intervals to facilitate comparison. The lowest rapidity interval, 
2 < y < 2.5, is not scaled. 

with the data is very good, the exception being the most forward rapidity measurement of _D*+ where the calculation 
is above the data. Interestingly, while the normalization of the D* and calculations are rather similar over all of 
the rapidity intervals, compare Figs.lSJb) and (c), there is a significant drop in the measured D* cross section at low 
PT between 3.5 < y < 4 and 4 < y < 4.5 that is not reproduced inthe calculations. 

We have shown that the calculated uncertainties on the total charm cross section can be considerably reduced by 
fitting the data with a next-to-leading order calculation. When the same fit parameters are used to calculate the 
leptons from heavy flavor decays in the FONLL approach, the results are still in agreement with the data. 

III. QUARKONIUM PRODUCTION IN THE COLOR EVAPORATION MODEL 

We now turn to a treatment of quarkonium production within this same framework. Perhaps the simplest ap- 
proach to quarkonium production is the color evaporation model (GEM) which treats heavy flavor and quarkonium 
production on an equal footing. The GEM was first discussed some time ago [13, and has enjoyed considerable 
phenomenological success when applied at next-to-leading order in the total cross section and leading order in the 
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quarkonium pT distribution [3^ [s^] ■ 

In the CEM, the quarkonium production cross section is some fraction, Fc, of all QQ pairs below the HH threshold 
where H is the lowest mass heavy-flavor hadron. Thus the CEM cross section is simply the QQ production cross 
section with a cut on the pair mass but without any constraints on the color or spin of the final state. The color 
of the octet QQ state is 'evaporated' through an unspecified process which does not change the momentum. The 
additional energy needed to produce heavy-flavored hadrons when the partonic center of mass energy, a/I, is less than 
2mH, the HH threshold energy, is nonperturbatively obtained from the color field in the interaction region. Thus 
the quarkonium yield may be only a small fraction of the total QQ cross section below 2mH- At leading order, the 
production cross section of quarkonium state C in a pp collision is 

o'c™(s««) =-^C' V / ds dxidx2 ff{xi,n%) f^{x2,^l%) ct,j(s, , (5) 

where ij = qq or gg and <Jij{s) is the ij — ^ QQ subprocess cross section. 

The fraction Fc must be universal so that, once it is fixed by data, the quarkonium production ratios should be 
constant as a function of -y/ii y a-nd pt- The actual value of Fc depends on the heavy quark mass, m, the scale 
parameters, the parton densities and the order of the calculation. It was shown in Ref. Q that the quarkonium 
production ratios were indeed relatively constant, as expected by the model. In addition, Ref. [s^ showed that the 
data on the J/ip and open charm cross sections as a function of ^/s in hadroproduction and W-yN in photoproduction 
have the same energy dependence. 

The data we use to obtain Fc for the J/ip are from the compilation by Maltoni et al. [s^. The data range from 
fixed-target experiments with center-of-mass energy 6.8 <_v/s < 41.6 GeV [sol-iSal to data from the CERN ISR at 
= 23 [13, 30 Hp, 30.6 [H], 31 JsS, 52 [13, [Si, 52.4 [11, 53 [l^JM, 62.7 [Ml, and 63 [13, [11 GeV. Data from 
the PHENIX experiment at RHIC 62] are also used. The ISR datalH^eil are all from pp measurements, as are the 
data from Refs. [13, Ek- Data from single nuclear targets include Bejliriii, 111, [Hi JLi jii^ H, [ii] , Si [Hi], Fe 
[49^ , Au [5^ , and Pt |43j| . Other experiments took data on multiple nuclear targets [H, I50l - l52l . [Hq . Both the total 
forward cross section (^f > 0) f39'- 56, and the cross section times the branching ratio to lepton pairs, Bu, at 
y ^ 0, Biida / dy\y=o (4 I |42|, |47|-|52„ ,54,-62] were reported. Several of the ISR experiments [57i - l59l. leij only provided 
the cross section at ?/ = 0, likely due to their limited phase space coverage. In cases where the total cross section 
was reported, the uncertainty provided was on the level of 40%. We note that several detectors have taken data at 
the same energy and with the same target but reported results with different experiment numbers that diverge by 
more than one standard deviation. For example, the p+C results reported by E331 [44| and E444 [4^ using the same 
apparatus, u = 256 ± 30 and 166 ± 23 nb/nucleon respectively, differ by more than two standard deviations. 

Maltoni and collaborators corrected prior measurements using up-to-date values of the J/ip branching ratios to 
IX^ IJL~ and e+e~ [63| a nd, when appropriate, averaged the results on multiple nuclear targets assuming UpA/cTpp = A"" 
with a — 0.96 ± 0.01 64| at ^ 0, obtained with an 800 GeV proton beam [s^]. The A dependence was assumed 
to be independent of center-of-mass energy. However, a recent reanalysis of these data, assuming a combination of 
shadowing and absorption effects on J production, found that, at 0:^ ^ 0, the absorption cross section decreases as 
a function of incident energy whether or not the data were corrected for shadowing effects j65j] . Later measurements 
withpiab — 158 GeV obtained an absorption cross section consistent with the predicted extrapolation [6^. In addition, 
effects arising from modifications of the parton densities in the nucleus that may be present in the data depend on 
the magnitude of A and have not been taken into account in the averaging. Thus the effective a, which includes all 
relevant nuclear effects, likely depends on incident energy. 

We have fit Fc to both the full data set as well as to more limited sets. Our final result is based on the total cross 
section data with only p, Be, Li, C, and Si targets respectively. In this way, we avoid uncertainties due to ignoring 
any cold nuclear matter effects which are on the order of a few percent in light targets. We also restricted ourselves 
to the forward cross sections only, rather than include the Buda / dy\y=o data in the fits. The rapidity distributions 
calculated in the MNR code are subject to fluctuations about the mean, even with high statistics calculations. The 
total cross sections, not subject to these fluctuations, are thus more accurate. 

Our calculations use the NLO QQ code of Mangano et al. (MNR) ^] with the 2m/f mass cut in Eq. ([5]), as 
described in Refs. [1, 111]. Because the NLO QQ code is an exclusive calculation, we take the mass cut on the 
invariant average over kinematic variables of the c and c. Thus, instead of deflning and ^.r relative to the quark 
mass, they are defined relative to the transverse mass, ijlf,r oc rriT = \J vrP- -\- pl^ where pt is that of the QQ pair. 



p\ = ^.h{p\^+p\_). 

We use the same values of the central charm quark mass and scale parameters as in the previous section, both for 
the fiducial parameter sets and for the best fit values, to obtain the normalization Fc- We fit Fc to the J/i/; data 
for both the fiducial FONLL sets (central value {m, I'm, ^n/'m) — (1.5 GeV, 1, 1)) and the fit results obtained in 
the previous section (central value {m, ^p/m, fiu/m) — (1.27 GeV, 2.1, 1.6)). We determine Fc only for the central 
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parameter set in each case and scale all the other calculations for that case by the same value of Fc to obtain the 
extent of the J/-0 uncertainty band employing Eqs. ([2]) and ([3]). 

We find Fc = 0.040377 for the central result, {m, fip/m, fin/m) = (1.5 GeV, 1, 1), and Fc = 0.020393 for the 
central CTIO result with {m, iip /m, ^n/m) = (1.27 GeV, 2.1, 1.6). A significantly larger value of Fq is necessary for 
the larger quark mass since the fraction of the total charm cross section remaining after the mass cut is smaller. The 
results for the energy dependence of the forward inclusive J/ip cross section are shown in Figs. |9] and [10] for the central 
mass values of 1.5 GeV and 1.27 GeV respectively. The uncertainty bands are shown on the left-hand sides of these 
figures while the individual parameter sets contributing to the bands are shown on the right-hand sides. 
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FIG. 9: (Color online) The forward J/t/i cross sections calculated with the FONLL fiducial parameter set with m = 1.5 
GeV. (a) The uncertainty band obtained employing the FONLL parameter set. The solid red curve is the central value 
while the limits of the uncertainty band are given by the dashed magenta curves. The dotted blue curve is a result with 
(m, hf /rriT, fiR/niT) = (1-2 GeV, 2, 2). (b) The solid red curve is the central value (m, ^p/niT, hr/itit) = (1.5 GeV, 1,1). The 
solid blue and magenta curves outline the mass uncertainty with (1.3 GeV, 1, 1) and (1.7 GeV, 1, 1) respectively. The dashed 
curves are associated with fj,/mT ~ 0.5: {nF/niT, fJ^R/rriT) = (1,0.5) blue; (0.5,1) magenta; and (0.5,0.5) red. The dotted 
curves are associated with fi/rriT ~ 2: {^f /mr , fJ.R/mT) = (1,2) blue; (2,1) magenta; and (2,2) red. 

The most obvious result in Fig. [HJa) is that there is no well defined lower limit on the total cross section with 
the fiducial parameter set, only an upper limit. The reason is apparent from Fig. [UJb): the combined differences 
in the minimum values of the masses and scales added in quadrature are larger than the central value for ^/s < 63 
GeV. When the fiducial parameter sets are applied to the CEM calculation of J/tp production, the upper limit of the 
charm quark mass, m = 1.7 GeV, gives a very narrow invariant mass interval for the CEM calculation in Eq. ([5]), 
from 2m — 3.4 GeV to 2mD = 3.86 GeV. The difference between the results with different quark masses is more 
pronounced at low center of mass energies while the energy dependence of the calculations with different values of 
m begins to converge at large y/s. Indeed, the 'hump' in the upper limit of the fiducial uncertainty band is due to 
the slower growth of {m, fip /niT, fJ-R/niT) = (1.3 GeV, 1,1) relative to (1.5 GeV, 1, 0.5) for ^/s > 400 GeV, note the 
crossing of the solid and dashed blue curves in Fig. [DJb). We also note that the fiducial set does not give very good 
agreement with the total J/tp cross section reported by CDF [6i| since the calculated ^/s dependence is too slow to 
match the measured growth of the forward cross section. 

The best fit band, shown in Fig. fTUTa). on the other hand, gives very good agreement with the J/4> data over the 
entire energy range, even for the CDF cross section, not included in the fit. The data are almost all encompassed by 
the width of the band. Now, as was the case for the total charm cross section, the uncertainty due to the quark mass 
dominates over that due to the scale choice for ^/s < 200 GeV. 

We now turn to the J/ip rapidity and pr distributions, shown in Fig. [11] for y/s = 200 GeV and Fig. [T^] for y/s — 7 
TeV. At leading order in the total cross section, the QQ pair px is zero. Thus, while our calculation is next-to-leading 
order in the total cross section, it is leading order in the quarkoniumpT distributions. In the exclusive NLO calculation 
[67j both the Q and Q variables are integrated to obtain the pair distributions, recall fip.R oc tot- 
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FIG. 10: (Color online) (a) The uncertainty band on the forward J/i/; cross section calculated based on the cc parameter fit 
in Fig. [2l^d). The dashed magenta curves and dot-dashed cyan curves show the extent of the corresponding uncertainty bands. 
The dashed curves outline the most extreme limits of the band, (b) The components of the uncertainty band. The central value 
{m, np/mT, IJ^R/mr) — (1.27 GeV, 2.10, 1.60) is given by the solid red curve. The solid blue and magenta curves outline the 
mass uncertainty with (1.18 GeV, 2.10, 1.60) and (1.36 GeV, 2.10, 1.60) respectively. The dashed curves outline the lower limits 
on the scale uncertainty: {fiF /rriT, j-iR/niT) = (2.10,1.48) blue; (1.25,1.60) magenta; and (1.25,1.48) red. The dotted curves 
outline the upper limits on the scale uncertainty: {^f /mr , i-iR/mT) ~ (2.10,1.71) blue; (4.65,1.60) magenta; and (4.65,1.71) 
red. The upper and lower dot-dashed cyan curves correspond to {jj.F/mT, fJ-R/mr) = (4.65, 1.48) and (1.25,1.71) respectively. 



Results on open heavy flavors indicate that some level of transverse momentum broadening is needed to obtain 
agreement with the low px data. This is often done by including some intrinsic transverse momentum, fc^, smearing to 
the initial-state parton densities. The implementation of intrinsic kx in the MNR code is not handled in the same way 
as calculations of other hard processes due to the nature of the code. In the MNR code, the cancellation of divergences 
is done numerically. Since adding additional numerical Monte-Carlo integrations would slow the simulation of events, 
in addition to requiring multiple runs with the same setup but different intrinsic kr kicks, the kick is added in the 
final, rather than the initial, state. In Eq. ([5]), the Gaussian function gp{kT), 

9p{kT) = -JTJT- cxp(-fc^/ (fc^)p) , (6) 

Tl[KrpJp 

[rij , multiplies the parton distribution functions for both hadrons, assuming the x and kx dependencies in the initial 
partons completely factorize. If factorization applies, it does not matter whether the kx dependence appears in the 
initial or final state if the kick is not too large, as described below. In Ref. [lH, {kx)p = 1 GeV^ was found to best 
describe fixed-target charm production. 

In the code, the QQ system is boosted to rest from its longitudinal center-of-mass frame. Intrinsic transverse 
momenta of the incoming partons, fcyi and fcT2j are chosen at random with k^^ and k^^ distributed according to 
Eq. (|ni). A second transverse boost out of the pair rest frame changes the initial transverse momentum of the QQ pair, 
Pt, to Pt + ^ti +kT2- The initial k^ of the partons could have alternatively been given to the entire final-state system, 
as is essentially done if applied in the initial state, instead of to the QQ pair. There is no difference if the calculation 
is leading order only but at NLO an additional light parton can also appear in the final state so the correspondence 
is not exact. In Ref. [7l|, the difference between the two implementations is claimed to be small if fc|, < 2 GeV^. We 
note that the rapidity distribution, integrated over all pT, is unaffected by the intrinsic kr- 

The effect of the intrinsic kx on the shape of the J/-0 px distribution can be expected to decrease as y/s increases 
because the average pr of the J/i/' also increases with energy. However, the value of (/cy) may increase with ^/s. We 
can check the energy dependence of {k'^) by the shape of the J/i/' pr distributions at central and forward rapidity 
at RHIC. We find that {k^) = 1 + (1/12) ln(yi/20) w 1.19 GeV^ at = 200 GeV agrees weU with the J/V Pt 
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FIG. 11: (Color online) The J/ip rapidity distribution (a) and the midrapidity (b) and forward rapidity (c) pt distributions and 
their uncertainties calculated with the same parameters as Fig. [Hd). The results are compared to PHENIX pp measurements 
at ^/s = 200 GeV [73|. The correlated and uncorrelated systematic errors in Ref. [T^] are added in quadrature. No additional 
scaling factor has been applied. The solid red curve shows the central value while the dashed magenta curves outline the 
uncertainty band. A (fey) kick of 1.19 GeV'^ is applied to the pr distributions, as discussed in the text. 

distributions measured by PHENIX at both midrapidity and forward rapidity, see Fig. 1111 The rapidity distributions, 
as well as the pr distributions in the two rapidity regions, all agree well with the J/ijj cross sections calculated with the 
central set of parameters. Only the low part of the forward rapidity pt distribution is somewhat underestimated. 
The integrated forward cross section is about 50% lower than the midrapidity value. In addition, the pt distribution 
falls off faster at high pT in the forward rapidity region. 




FIG. 12: (Color online) The J/ij) rapidity distribution (a) and the midrapidity, \y\ < 0.9 (b), and forward rapidity, 2.5 < y < 4 
(c) Pt distributions at ^/s = 7 TeV and their uncertainties calculated with the same parameters as in Fig. JHd). The results 
are compared to the ALICE rapidity distribution as well as the mid and forward rapidity pt distributions [73|. No additional 
scaling factor has been applied. The solid red curve shows the central value while the dashed magenta curves outline the 
uncertainty band. A (fcy) kick of 1.49 GeV^ is applied to the pt distributions, as discussed in the text. 

The ALICE 7 TeV px distributions, shown in Fig.H^b) and (c), include the ALICE rapidity cuts for the central and 
forward rapidity regions, \y\ < 0.9 and 2.5 < y < 4, respectively. The rapidity distribution at y/s = 7 TeV, Fig. [T^ a). 
is flat over several units of rapidity. Thus the integrated cross sections in the two rapidity intervals, normalized per 
unit of rapidity, are very similar. However, the forward rapidity pT distribution is still a stronger function of pt than 
the midrapidity distribution. 

Finally, the inclusive J/ip rapidity distribution and forward rapidity distribution at ^/s — 2.76 TeV are compared 
to the ALICE data in Fig. [131 Here the calculated rapidity distribution is not as broad and the agreement with the 
data is rather good although the midrapidity point remains high relative to the central value of the calculation. The 
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FIG. 13: (Color online) The J/ip rapidity distribution (a) and the forward rapidity, 2.5 < y < 4 pr distribution (b) at 
-^s = 2.76 TeV and their uncertainties calculated with the same parameters as in Fig. [lid). The results are compared to the 
ALICE rapidity distribution as well as the forward rapidity pr distribution [t^. No additional scaling factor has been applied. 
The solid red curve shows the central value while the dashed magenta curves outline the uncertainty band. A (fcy) kick of 1.41 
GeV^ is applied to the pr distributions, as discussed in the text. 

agreement of the calculated px distribution with the forward rapidity data is quite good with the exception of the 
lowest pt point where the calculated distribution turns over more quickly than the data. 

IV. SUMMARY 

We have narrowed the uncertainty band on the open heavy flavor cross section and, in so doing, have also provided 
a realistic uncertainty band on J/tjj production in the color evaporation model. The central result, m — 1.27 GeV, 
fip/m = 2.1 and /ifl/m = 1.6, is quite compatible with previous calculations using a 'by-eye' fit to the data with 
TO = 1.2 GeV, fip/m = fin/m = 2%^. 

While the fits have been made by comparing the calculated NLO charm production cross section to available data 
at fixed-target energies and at RHIC, they are in good agreement with the extracted total charm cross sections at the 
LHC. The same parameter set also provides good agreement with the distributions of single leptons from semileptonic 
heavy flavor decays at RHIC and the LHC. The limit on the width of the uncertainty band is now set by the uncertainty 
due to bottom quark production and decay. 

We have used the same fit parameters in the calculation of J/tp production in the color evaporation model and have 
thus provided the first uncertainty band on J/tp production in this approach. The energy dependence of the total 
J/ip cross section that results is a good match to the data up to collider energies. The pt distributions are also in 
good agreement with the data from RHIC and the LHC. In future work, we will use our new parameter set to place 
limits on the contribution of B meson decays to J/tA production and will also study cold nuclear matter effects on 
J/ip production. 
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